package scu.maqiang.eigenvalue;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T32;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class ElasticityEigenValue2DT32 {

	public static void main(String[] args) {
		
		String meshFile = "mesh2D\\SteelBracket.msh";
		Mesh2T3 mesh = new Mesh2T3();
		mesh.readMesh(MeshFormat.msh, meshFile);
		FES2T32 fs = new FES2T32(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		double E = 1.0, Nu = 0.3, shift = 2.0;
		
		
		fs.assembleStiff(new double[] {E, Nu}, BVPType.PLANESTRESS, A);
//		A.sort();
//		System.out.println(A);
//		System.exit(0);
		fs.assembleMass(new double[] {1.0}, BVPType.CONSIST_MASS, M);
		
		int N = 10;
		double[] x = new double[fs.GetNdof()];
		Arrays.fill(x, 1.0);
		double[][] disp = new double[2][mesh.getNv()];
		double[] D = new double[N];
		double[][] V = new double[N][fs.GetNdof()];
		String outputFilePreFix = "ElasticEigenValue";
		EigenSSolver solver = new EigenSSolver();
		solver.inversePowerTran_PLU(A, M, fs.getnBand(), x, shift, N, D, V);
		for(int i = 0; i < N; i++) {
			fs.extractComponent(V[i], disp);
			mesh.toTecplot(outputFilePreFix + (i+1) + ".dat", disp);
		}
		System.out.println(MVO.toString(D));
	}
}
